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Abstract 

We present a new algorithm for solving the real roots of a bivariate polynomial system E = [.fix, y), g{x, y)] with a finite 
number of solutions by using a zero-matching method. The method is based on a lower bound for bivariate polynomial 
system when the system is non-zero. Moreover, the multiplicities of the roots of E = can be obtained by a given 
neighborhood. From this approach, the parallelization of the method arises naturally. By using a multidimensional 
matching method this principle can be generalized to the multivariate equation systems. 
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1. Introduction 

Considering the following system: 



^^{.f{x,y),g{x,y)}, (1) 



we assume that f{x,y),g(x,y) e where Q is the field of rational numbers. We call the zero-dimension if the 

bivariate polynomial system ([T]) has a finite number of solutions. 

Real solving bivariate polynomial system in a real field is an active area of research. It is equivalent to finding 
the intersections of f{x,y) and g{x,y) in the real plane. The problem is closely related to computing the topology 
of a plane real algebraic curve and other important operations in non-linear computational geometry and Computer- 
Aided Geometric DesignlU fTslflsiflolflsil . Another field of application is the quantifier elimination 1 7i [l7tL_Jhere 
are several algorithms that tackle this problem such as the Grober basis method[ll9i the resultant method |23Q, the 
characteristic set method and the subdivision methodfl 21 1. However, the procedure of these techniques is very 



complicated. In this paper, we propose an efficient approach to remedy these drawbacks. 

In this paper, we propose a zero-matching method to solve the real roots of an equation system like ([TJ. The basic 
idea of zero-matching method is as follows: First projecting the roots of S to the x-axis, gives the roots {x\, - ■ ■ , x,,}, 
and the y-axis, gives the roots {yi, ■ • ■ ,yv}, respectively. Subsequently, for every root x,, and for every yj is back- 
substituted in f{x,y) and g(x,y). To that end, for some root x, there is the corresponding one or more roots yj to be 
determined satisfying S. The main contribution of our method is that how to determine the real roots of E = and 
the multiplicities of the roots. Moreover, our approach that has given solutions to this situation can be the design of 
parallelized algorithms. 
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In i^, Diochnos et al presented three algorithms to real solving bivariate systems and analyzed their asymptotic 
bit complexities. Among the three algorithms, the difference is the way they match solutions. The method of special- 
ized Rational Univariate Representation (rur) based on fast gcd computations of polynomials with coefficients in an 
extension field to achieve efficiency (hence the name g_rur) has the lowest complexity and performs best in numerous 
experiments. The grur method projects the roots to the x-axis and y-axis, for each x-coordinate a computes the 
GCD h(a, y) of the square-free parts of f(a, y) and g(a, y), and isolates the roots of h(a, y) = based on computations 
of algebraic numbers and the rur techniques. Our algorithm only uses resultant computation and real solving for 
univariate polynomial equations with rational coefficients. 

The hybrid method proposed by Hong et al[ 16] that projects the roots of E to the jc-axis and y-axis respectively and 
uses the improved slope-based Hansen-Sengupta to determine whether the boxes formed by the projection intervals 
contain a root of S. The numerical method only works for simple roots of S. When the system has multiple roots, the 
rur technique is used to isolate the roots. Compared with this method, our approach also computes two resultants of 
the same total degrees. However, our method is a complete one, their numerical iteration method needs to use the rur 
technique to find multiple roots. 

In||2j], Bekker et al presented a Combinatorial Optimization Root Selection method (hence the name cors) to 
match the roots of a system of polynomial equations. However, the method is only suitable for solving a small 
system of polynomial equations, and does not work for the multiple roots. Recently, Cheng et al\4] proposed a local 
generic position method to solve the bivariate polynomial equation system. The method can be used to represent the 
roots of a bivariate equation system as the linear combination of the roots of two univariate equations. Moreover, 
the multiplicities of the roots of the bivariate polynomial system are also derived. However, the method is very 
complicated to extend to solve the multivariate equation systems. Our method can solve the larger systems and easily 
generalize to the multivariate equation systems. 

The rest of this paper is organized as follows. In Section 2, we give some notations, a lower bound for bivariate 
polynomial equation if it is non-zero, and how to determine the root multiplicity. In Section 3, we propose the 
algorithm to real solving the bivariate polynomial system and give a detailed example. In section 4, we present some 
comparisons of our algorithm. The final section concludes this paper. 



2. Notations and main results 

2.1. Notations 

In what follows D is a ring, F is a commutative field of characteristic zero and F its algebraic closure. Typically 
D = Z, F = Q and F = Q. 

In this paper, we consider the zero-dimensional bivariate polynomial system as follows: 

0<i<« 0< i<m 
0<i<pO<j<q 

Throughout this paper, note that deg^ = max(n,p), degy = max(m,<7), = max(||/||i, \\g\\\), where the ||/||i and 
WgWi are the one norm of the vector (ooo, aoi, • ■ ■ , aom, ■■■ , a«o, ■ ■ • , a„„) and (bau boi,-- - , bog, ■■■ , bj^, ■■■ , bpq), so 
ll/lli - and \\g\\\ - E,Ey|/7,y|, respectively. M - max(||f||i, where the t{x) and T{y) are the no extraneous 

factors in resultant polynomial of Z. \L\ denotes that the bivariate polynomial system S has been assigned values to 
two variables. 

Let n be the projection map from the Z to the x-axis: 

TT : ^ R, such that n(x,y) = x. (3) 

For a zero-dimensional system 2 defined in (|2]l, let f(x) e Q[x] be the resultant of f{x,y) and g{x,y) with respect to y: 

t{x)^RQSy{f{x,y),g{x,y)). (4) 
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Since S is zero-dimensional, we have t(x) ^ 0. Then :7r(V(E)) c V(f(x)), where V(/i , ■ ■ ■ , /„,) is the set of common 
real zeros of = 0. If t(x) is irreducible, then denote the highest degree by deg,. Let the real roots of t{x) = be 

ai < Q'2 < • • • < «„. (5) 

By using the same method, let T{y) € Q[y] be the resultant of f{x, y) and g{x, y) with respect to x: 

T(y)^ Res Af(x,y),8(x,y)). (6) 

If T(y) is irreducible, then denote the highest degree by degr- Let the real roots of T(y) be as follows: 

/3i </32<--- </3,. (7) 

We observe that the above projection map may generate extraneous roots. Fortunately, we can easily discard these 
extraneous factors by computing the determinant of the sub-matrix of the coefficient matrix. Moreover, if the resultant 
is irreducible, then it is no extraneous factors. However, when the resultant is reducible, it may suffer from the 
extraneous factors. The method of removing extraneous factors mentioned can be adapted to the resultant for the 
bivariate polynomial system |27]. It is the following theorem to remove the extraneous roots. 

Theorem 2.1. 2 is defined as in (|2|. If the resultants of ^ for one variable is reducible, denoted by tern, then the 
resultant of bivariate polynomial system is the only some irreducible factors in which the other variable appear. 

Proof. The proof can be given similarly to that in Proposition 4.6 of Chapter 3 of |8]. □ 

2.2. A lower bound for if 1, ^ 

The purpose of this subsection is to prove the following theorem. 

Theorem 2.2. Z is defined as in Let a,f3 be two approximate real algebraic numbers. Denote by the integer 
s — deg, ■ degT, and N as above. If |S| ^ 0, then 

|I| > N^-'M-'--\ (8) 

where c is the constant satisfying certain conditions, |S| is the following two cases: 

(a) Iff{a,l5) = org{a,p) = 0, then = max{|/(a,/3)|, \g{a,p)\]; 

(b) Iff(a,f5) + andg(a,fi) + 0, then \Y\ = min{|/(a,/?)|, |g(a,yS)|}. 

Before giving the proof of theorem lZ2l we recall two lemmas: 

Lemma 2.1. ( ^Idi . lemma 3) Let ay, . . . ,aq be algebraic numbers of exact degree of d\, . . . ,dq respectively. Define 
D — [Q(q;i, . . . , Uq) : Q]. Let P e Z[xi, . . . ,Xq] have degree at most Nh in Xh(l < h < q). IfP(ai, . . . , a^) + 0, then 

q 

\P{a,,. . . , aq)\ >\\ P |||-« Yl Miahr""'-''", 

h=l 

where the M(ah) is the Mahler measure of a^. 

Proof See the Lemma 4 of |20]. □ 

Lemma 2.2. Let a be an algebraic number Denote by the M{a) of the Mahler measure of a. If P is a polynomial 
over Z, then 

M(a) < \\P\\i. 

Proof For any polynomial P - 'Z,'!=o Pi ^ of degree d with the all roots cr*'', • ■ • , cr'-''\ we define the measure 
M{P) by 

M(P) = |pdintimax{l,|cr«|). 

The Mahler measure of an algebraic number is defined to be the Mahler measure of its minimal polynomial over Q. 
We know from Landau (" lll4ll . p. 154, Thm. 6. 31) that for each algebraic number a 

M(a) < \\P\\2, 

where ||f||2 = (lif=o Ip/P)'^^- It is very easy to get that ||f II2 < lli^lli. This completes the proof of the lemma. □ 

3 



Now we turn to give the proof of Theorem l2.2l 



Proof. From the assumption of the theorem, since E is defined as in (|2]l. Let the pair {a, /3) be corresponding value to 
the variable x and y for S respectively. We have the following equations: 

f(a,/^) = Z ""'-j"'^' ^^^^ 

0</<n 0< j<m 
0<i<p 0< j<q 

At first, we consider the lower bound for the equation ( |9a] i. Define k - [Q(a,/3) : Q]. Denote by \f\=\f(a,/3)\, and r, f 
by the exact degree of algebraic numbers a, j3 respectively. From Lemma ( 12.11 ). if |/| ^ 0, then 

I/I > 11/11 j"*M(a)-*"/''M(/3)-'^"'/'. 

We observe that M{a) and M(j6) derive from t(x) and ^(y) respectively. From Lemma ( 12.2b . we can get the following 
inequality: 

M(a)<||f||i,M(/?)<||r||i. 

So we can obtain that 

I/I > ii/ir^^iifiir^'^niT'ii:""^'. (10) 

By using the same technique as above, we can obtain the lower bound for the equation (l9b] i. Denote by |g|=|g(a,y6)|. 
If \g\ + 0, then 

\g\>\\gr%t\v,'"^%w,'"'>\ (11) 

Since we have the following two cases: 

(a) \ff(a,fi) = org(a,/?) = 0, then |E| = max{|/(a,/?)|, \g(aM\\ 

(b) If /(a,^) + and g(,a,p) + 0, then |S| = min{|/(a,^)|, |^(a,yS)|}. 

Hence we are able to obtain the lower bound for the bivariate polynomial system. From the above assumption, we can 
get the following parameters: 

k = [Q(a,yS) : Q] < deg{t{x)\ ■ deg{T(y)] = deg, ■ degr, (12a) 
N = max{||/||i, \\g\\i],M^ max{||f||i, ||r||i), r = deg„t = degj. (12b) 

Combined with the equation ( I12al l and ( I12bl l. it is obvious that s - k and the constant c - ^ + + 1 . Finally, note 
that the constant c satisfies both cases. This proves the theorem. □ 



As the corollary of Theorem l2.2l we have 

CoroUay 2.1. Under the same condition ofTheorem \2.2\ if |!E| < A^'^'M"'^ then |E| - 0. We say that a is associated 
with Pfor the real root of 2 . Denote by the s — N^^^ M^'^'" for the rest of this paper 

Proof. The proof is very easy by contradiction. □ 
2.3. Root multiplicity 

The results of this subsection can be provided for the root multiplicity of S. We follow the approach and terminol- 
ogy of isll andlflUl. 

Let C/, Cg be /, g corresponding affine algebraic plane curves, defined by the equations S. Let / =< /, g > be the 
ideal that they generate in F[jt:,y], and so the associated quotient ring is J?l = F[jic,3']//. Let the distinct intersection 
points, which are the distinct roots of (S), be C/ n Cg c [S ij - (ai,[5j)}\<i<u,i<j<v 

The multiplicity of a point S ij is 

mult(S ij '. Cf C\ Cg) — dimf^Sij < °°> 
4 



where ^Sij is the local ring obtained by localizing J?l at the maximal ideal I -< x - ai,y - /3j >. 

If ^Sij is a finite dimensional vector space, then 5, y = {o:i,/3j) is an isolated zero of / and its multiplicity is called 
the intersection number of the two curves. The finite J?l can be decomposed as a direct sum J?l = yis,, ^ © " ' ' 
and thus diirifjl = 2"=i mult(Sij : C/ n Q). 
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Proposition 2.1. ( hllxl . Proposition 1 ) Let f,g& F[x, y] be two coprime curves, and let p €¥ be a point. Then 

mult(p : fg) > multip : f)mult{p : g), 
where equality holds if and only ifCf andCg have no common tangents at p. 

Proposition 2.2. Let us obtain the real roots of ^ — in Q and If the two matching pairs (ai,/3j) and 
(ai+i,Pj+i){for 1 < / < M, 1 < 7 < v) are satisfying S = 0, Ic; - q;,+i| < e and \/3j - /3j+i\ < s, then the {a-^Pj) 
is multiple root of 1, — 0. 



Proof. From Theorem l2.2l and Corollarv l2.1l it is obvious that E = if and only if 

HI < s. 

Therefore, the error controlling is less than s in numerical computation. Under the assumption of the proposition, we 
get Id-, - £!',■+ 1 1 < s and \/3j - /3j+i \ < s. So we are able to obtain that |q',- - q;,+i| = and -y6;+i| - int the truncated 
error This proves the proposition. □ 

From Corollarv l2.1l the two-tuple {a,/3) is the real root of S = 0. This method is called a zero-matching method. 
The technique is a posteriori method to match the solutions for the bivariate system. It can be generalized easily to 
real solving the multivariate polynomial systems. 



3. Derivation of the Algorithm 

The aim of this section is to describe an algorithm for real solving bivariate polynomial equations by using zero- 
matching method. We first find the parameters A^, c and s, then obtain the no extraneous factors t(x) and T(y) with the 
resultant elimination methods, and real solving two univariate polynomials, and finally match the real roots for the 
systems. 

5.7. Description of algorithm 

Algorithm 1 is to discard the extraneous factors from the resultant method, algorithm 2 is to obtain the solutions 
of bivariate polynomial systems. 

Algorithm 1 NoExtrRes(S, var) 

Input: {f(x,y),g{x,y)], var is one variable. 

Output: No extraneous factors resultant of S. 

1: tem <^ Res^ar{f(x,y),g(x,y)]; 

2: if tem is irreducible then 

3: return tem ; 

4: else 

5: tem <— Res ■ extraneous facotrs; 
6: return Res. 
1: end if 



Now we can give the algorithm|2]to compute the real roots for E = 0. 

The parallelization of the algorithm that we have just described can be easily done because it performs the same 
computations on different steps of data without the necessity of communication between the processors. Observe that 
the Step 1 and Step 2, Step 6 and Step 7 of the algorithm can be easily paralleled, respectively. 

Now we get a theorem about the computational complexity of the whole algorithm. 
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Algorithm 2 zmm(S) 

Input: E = {f{x,y), g(x,y)] is a zero-dimensional bivariate polynomial system. 

Output: A set for the real roots of E = 0. 
1: Project on the x-axis such that t{x) = Resy(f(x,y),g{x,y)); 
2: Project on the y-axis such that T{y) = Resx{f{x,y), g(x,y)y, 
3: Discard the extraneous factors from t(x) and T(y) by using Algorithm[T] 
4: Find the parameters and s, and Compute c according to the Theorem l2.2l 
5: Obtain the lower bound s by Corollarv l2.1l 

6: Solve the real roots of the resultant t(x) for the set Sx = {ai , a2, ■ ■ ■ ,a„]; 
1: Solve the real roots of the resultant T{y) for the set Sy - {y6i ,j02, ■ • ■ ,Pv]\ 

8: Match the real root pair to get the solving set S = {(ai,/3j), 1 < / < m, 1 < j < v} by Corollarv l2.1l 
9: Check the root multiplicity of the set S by Proposition |22] 



Theorem 3.1. Algorithm\2\works correctly as specified and its complexity includes as follows: 

(a) 0(dT + dlgd) for computation of real solving univariate polynomial, where d is the degree of corresponding 
polynomial, t = 1 + max,<^ lg\ai\ and a,- is the coefficients. 

(b) 0(uv)for matching the solutions of bivariate polynomial system. 

Proof. Correctness of the algorithm follows from theorem l272l 

(a) The number of arithmetic operations required to isolate all real roots is the number of real root isolation of uni- 
variate polynomial by using subdivision-based Descartes' rule of sign. Using exactly the same arguments we know 
that they perform the same number of steps, that is 0(dT + dlgd). 

(b) As indicated before, the problem of matching the real roots of polynomial system mainly relies on the scale of 
solutions of every variable, respectively. □ 

3.2. A small example in detail 

Example 3.1. We propose a simple example f(x,y) — x^-y^-3 andg(x,y) — 3x^—2y^ — l to illustrate our algorithms. 
Step 1: tix) = 4 * - 45 * jc"* -H 114 * - 109; 
Step 2: T(y) = (-2 * / -H 8 -I- 3 * y^; 

Step 3: Discard the extraneous factors T{y) = -3*}'^ — 8h-2* y^; 
Step 4: Obtain the parameters N — 5, c — 2, s — 4; 
Step 5: Obtain the lower bound s = .1280 X lO"**; 

Step 6: Solve the real roots of the resultant t(x) for the set Sx - {-2.858288520, 2.858288520); 
Step 7: Solve the real roots of the resultant T(y)for the set Sy — {2.273722337); 

Step 8: Combine the the pairs from Sx and Sy respectively. Substitute the pairs into ^for variables x andy, 
determine whether less than the lower bound s, finally we find that the pairs S — {{x — -2.858288520,}' = 
2.273722337), {x = 2.858288520, y = 2.273722337)) are the solutions for S; 

Step 9: The multiplicity of the root of the system is one. 

3.3. Generalization and applications 

As for the generalization of the algorithm to real solving the multivariate equation systems case, we have to say that 
the situation is completely analogous to the bivariate case. However, its key technique is to transform the multivariate 
polynomial equations to the corresponding univariate polynomial equations. We can consider the Dixon Resultant 
Method to break this problem [^. However, we observe that how to improve the projection algorithm in resultant 
methods is the significant challenge. 
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Moreover, our algorithm is applicable for rapidly computing the minimum distance between two objects collision 
detection jzsll . This also enables us to improve the complexity of computing the topology of a real plane algebraic 
curve |9|]. 

4. Some comparisons 

We have implemented the above algorithms as a software package zmm in Maple 12. For problems of small size 
Uke the example of Section 3, any method can obtain the solutions in little time. But when the size of the problems is 
not small the differences appear clearly. Extensive experiments with this package show that this approach is efficient 
and stable, especially for larger and more complex bivariate polynomial systems. 

We compare our method with lgp |4], Isolate i23ll . discoverer |24], and grur lgp is a software package for 
root isolation of bivariate polynomial systems with local generic position method. Isolate is a tool to solve general 
equation systems based on the Realsolving c library by Rouillier discoverer is a tool for solving semi-algebraic 
systems, grur is a tool to solve bivariate equation systems. The following examples run in the same platform of 
Maple 12 under Windows and amd Athlon(tm) 2.70 ghz, 2.00 gb of main memory. We did three sets of experiments. 
The precision in these experiments is set to be high. In three tables, where ' ?' represents that the computation is not 
finished. 

In Table 1 the results are given both / and g are randomly generated dense polynomials with the same degree and 
with integer coefficients between -20 and 20. The command of Maple is as follows; 
randpoly([x,y],coeffs — rand(— 20. .20), dense, degree — 10). 



Table 1 : time for computing dense bivariate polynomials with no multiple roots 



system 


deg 


solutions 


Average Time(sec) 


/ 


8 


zmm 


LGP 


Isolate 


DISCOVERER 


grur 


SI 


4 


7 


2 


0.031 


0.031 


0.047 


0.313 


2.734 


S2 


6 


8 


6 


0.415 


1.328 


0.500 


1.828 


247.203 


S3 


7 


8 


6 


1.204 


2.734 


1.500 


7.047 


382.640 


S4 


8 


9 


6 


4211 


8.906 


4.672 


20.437 


2714.438 


S5 


9 


10 


2 


4070 


8.485 


4.687 


89.235 


1645.312 


S6 


10 


7 


6 


1.805 


3.860 


2.109 


22.250 


978.421 


S7 


10 


11 


4 


21.078 


43.734 


22.828 


7 


7 


S8 


12 


11 


2 


26.945 


54.969 


29.094 


7 


7 


S9 


12 


13 


4 


118.266 


241.734 


123.469 


7 


7 


SIO 


13 


11 


1 


15.446 


31.485 


17.796 


7 


7 


Sll 


14 


10 


8 


63.914 


200.828 


68.594 


7 


7 



In Table 2 the results are given both / and g are randomly generated sparse polynomials in the same degree, with 
sparsity default, and with integer coefficients between -20 and 20. The command of Maple is as follows: 
randpoly([x,y],coeffs — rand{— 20.. 20), sparse, degree — 10). 
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Table 2: time for computing sparse bivariate polynomials with no multiple roots 



system 


deg 


solutions 


Average Time(sec) 


/ 


8 


ZMM 


LGP 


Isolate 


DISCOVERER 


GRUR 


SI 


5 


6 


1 


0.015 


0.032 


0.015 


0.141 


1.032 


S2 


6 


7 


3 


0.040 


0.062 


0.047 


0.188 


5.375 


S3 


7 


5 


3 


0.024 


0.047 


0.047 


0.265 


2.688 


S4 


8 


6 


5 


0.031 


0.031 


0.047 


0.094 


1.031 


S5 


9 


8 


2 


0.047 


0.172 


0.078 


1.828 


51.000 


S6 


10 


11 


3 


0.063 


0.297 


0.125 


0.656 


11.110 


S7 


11 


9 


2 


0.164 


0.609 


0.375 


3.938 


877.875 


S8 


12 


13 


2 


1.141 


2.593 


1.453 


6.703 


1607.719 


S9 


13 


11 


4 


2.508 


5.344 


2.969 


7 


7 


SIO 


15 


17 


1 


0.532 


1.234 


1.266 


7 


7 


Sll 


20 


17 


4 


18.180 


39.688 


20.235 


7 


7 



In Table 3 the results are given is done with polynomial systems with multiple roots. We randomly generate a 
polynomial h{x,y,z) and take fix,y) = Res~{h, h.), g(x,y) - fy{x,y). Since f(x,y) is the projection of a space curve 
to the xy-plane, it most probably has singular points and f = g - Q is an equation system with multiple roots. The 
command of Maple is as follows: 

h :— mndpoly([x,y,z],coeffs — rand{-5 ..5), degree — 5);/ := re sultant{h, dif f (h, z), z)', g '■— diff{f,y). 



Table 3: time for computing bivariate polynomials with multiple roots 



system 


deg 


solutions 


Average Time(sec) 


/ 


g 


ZMM 


LGP 


Isolate 


DISCOVERER 


GRUR 


SI 


3 


2 


2 


0.016 


0.016 


0.016 


0.016 


0.062 


S2 


4 


3 


2 


0. 


0.032 


0.031 


0.016 


0.094 


S3 


4 


6 


7 


0.024 


0.016 


0.047 


0.109 


1.109 


S4 


5 


4 


3 


0. 


0.016 


0. 


0.016 


0.109 


S5 


6 


5 


2 


0.015 


0. 


0. 


0.016 


0.063 


S6 


9 


8 


2 


0.016 


0.046 


0.032 


0.015 


0.063 


S7 


12 


11 


3 


0.109 


0.234 


0.187 


0.063 


0.094 


S8 


13 


12 


7 


2.875 


137.641 


3.141 


1.328 


207.094 


S9 


14 


13 


4 


0.860 


2.891 


0.953 


0.141 


0.3110 


SIO 


19 


18 


1 


0.672 


1.547 


0.797 


22.156 


1520.812 


Sll 


16 


15 


5 


7.945 


27.047 


9.000 


7 


7 



From the Table 1, 2 and 3, we have the following observations. 

In the first two cases, the equations are randomly generated and hence may have no multiple roots. For systems 
without multiple roots, zmm is the fastest method, which is significantly faster than lgp and Isolate. Both zmm and lgp 
compute two resultants and isolate their real roots, lgp is slow, because the polynomials obtained by the shear map are 
usually dense and with large coefficients discoverer and grur generally work for equation systems with degrees 
not higher than ten within reasonable time. 

For systems with multiple roots, in the sparse and low degree cases, all methods are fast. Note that our method 
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is quite stable for equation systems with and without muhiple roots, lgp and Isolate are also quite stable, but slower 
than ZMM for bivariate equation systems. 

We also observe that all methods spend more time with sparse and dense polynomials than polynomials with 
multiple roots in the same high degree. This phenomenon needs further exploration. 

Remark 4.1. Of course, we should mention that discoverer and Isolate can be used to solve general polynomial 
equations and even inequalities. Here our comparison is limited to the bivariate case. In further work, we would like 
to consider solving multivariate polynomial equations. 

Remark 4.2. As is well known, the parallel algorithm is well suited for the implementation on parallel computers that 
allows the increase of the calculation speed. If our algorithm have been fully parallelized by using a large enough 
number of processors for each case, the real solutions of all the examples will have been computed in a couple of 
seconds. 

5. Conclusion 

In this paper, we propose a zero-matching method to real solving bivariate polynomial equation systems. The 
basic idea of this method is to find the lower bound for bivariate polynomial system when the system is non-zero. 
Moreover, we provide an algorithm for discarding extraneous factors with resultant and show how to construct a 
parallelized algorithm for real solving the bivariate polynomial system. An efficient method for multiplicities of the 
roots is also derived. The complexity of our method has increased steadily with the growth of bivariate polynomial 
system. Extensive experiments show that our approach is efficient and stable. The result of this paper can be extended 
to real solving of bivariate polynomial equations with more than two polynomials by using the resultant method. 
Furthermore, our method can be generalized easily to multivariate polynomial systems. 
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